####################################################################
##                                                                   
##      PROJECT TITLE:    Drummond et al. (2018) replication               
##      PROJECT AUTHOR:   SAM S. ROWAN                        
##      EMAIL:            SAM.ROWAN@NUFFIELD.OX.AC.UK         
##                                                                  
##      DESCRIPTION:      R code to make figures in Pitfalls in comparing Paris pledges
##                                                                  
######################################################################

library(haven)
library(tidyverse)
library(stargazer)

# rm(list=ls())
# Set working directory to the folder containing (code, data, figures, table)
setwd("/Users/sr/Dropbox/Projects/DHSP/dhsp-july")

#####################################################################

#### Figure 1: GHG trajectories and targets ####

ghg <- read_dta("data/primap-1988panel.dta")

## Canada
ghg %>% 
  filter(iso == "CAN") %>% 
  mutate(ghg = primap_l_ghg/1000) -> canada

ggplot(canada, aes(year, ghg)) + 
  # Add historical emissions
  geom_line(size=1.75) +
  # Add lines to signify the targets
  geom_segment(x = 1990, xend=2030.5, 
               y = canada$ghg[canada$year==1990], yend = canada$ghg[canada$year==1990],
               col="blue", linetype=2, size=1.5) +
  geom_segment(x = 1990, xend=2029.5, 
               y = canada$ghg[canada$year==2005], yend = canada$ghg[canada$year==2005],
               col="red", linetype=2, size=1.5) +
  geom_segment(x=2030.5, xend=2030.5, 
               y=canada$ghg[canada$year==1990], yend=canada$ghg[canada$year==2005]*.7,
               col="blue", linetype=2, size=1.5) +
  geom_segment(x=2029.5, xend=2029.5, 
               y=canada$ghg[canada$year==2005], yend=canada$ghg[canada$year==2005]*.7,
               col="red", linetype=2, size=1.5) +
  # Add INDC marker
  geom_point(aes(x=2030, y=canada$ghg[canada$year==2005]*.7), pch=4, col="red", size=4) +
  # Add text labels
  annotate("text", x = 2029, y = canada$ghg[canada$year==2005]*.7, 
           label="INDC target", adj=1, col="red", size=6) +
  annotate("text", x = 2029, y = canada$ghg[canada$year==2005], 
           label = "30% reduction\nfrom 2005 levels", adj=1, size=6) +
  annotate("text", x = 2029, y = canada$ghg[canada$year==1990], 
           label = "17% reduction\nfrom 1990 levels", adj=1, size=6) +
  # Settings
  coord_cartesian(xlim=c(1990,2031), ylim = c(500,950)) +
  labs(x="Year", y="GHG emissions (MtCO2)") +
  theme(legend.position = "none",
        panel.background = element_rect(fill=NA), 
        axis.text = element_text(size = rel(1.25)),
        panel.grid.major.y = element_line(colour = "grey90"),
        panel.ontop = FALSE,
        text=element_text(size=18),
        axis.line = element_line(size = .5, colour = "black"))
ggsave("figures/canada-traj.pdf", width=8, height = 6, units = "in")

# Russia
############

ghg %>% 
  filter(iso == "RUS") %>% 
  mutate(ghg = primap_l_ghg/1000) -> russia

ggplot(russia, aes(year, ghg)) + 
  # Add historical emissions
  geom_line(size=1.75) +
  # Add lines to signify targets
  geom_segment(x = 1990, xend=2030.5, y = russia$ghg[russia$year==1990], yend = russia$ghg[russia$year==1990],
               col="red", linetype=2, size=1.5) +
  geom_segment(x = 1990, xend=2029.5, y = russia$ghg[russia$year==2005], yend = russia$ghg[russia$year==2005],
               col="blue", linetype=2, size=1.5) +
  geom_segment(x=2030.5, xend=2030.5, y=russia$ghg[russia$year==1990], yend=russia$ghg[russia$year==1990]*.725,
               col="red", linetype=2, size=1.5) +
  geom_segment(x=2029.5, xend=2029.5, y=russia$ghg[russia$year==2005], yend=russia$ghg[russia$year==1990]*.725,
               col="blue", linetype=2, size=1.5) +
  # Add marker points for the target range, and connect them
  geom_point(aes(x=2030, y=russia$ghg[russia$year==1990]*.7), pch=4, col="red", size=4) +
  geom_point(aes(x=2030, y=russia$ghg[russia$year==1990]*.75), pch=4, col="red", size=4) +
  geom_segment(x=2030, xend=2030, y=russia$ghg[russia$year==1990]*.71, yend=russia$ghg[russia$year==1990]*.74,
               col="red", linetype=1, size=1) +
  # Add text labels
  annotate("text", x = 2029, y = russia$ghg[russia$year==1990]*.725, 
           label="INDC target range", adj=1, col="red", size=6) +
  annotate("text", x = 2029, y = russia$ghg[russia$year==1990], 
           label = "25-30% reduction \nfrom 1990 levels", adj=1, size=6) +
  annotate("text", x = 2029, y = russia$ghg[russia$year==2005], 
           label = "9-17% increase \nfrom 2005 levels", adj=1, size=6) +
  # Settings
  coord_cartesian(xlim=c(1990,2031), ylim = c(2000,3750)) +
  labs(x="Year", y="GHG emissions (MtCO2)") +
  theme(legend.position = "none",
        panel.background = element_rect(fill=NA), 
        axis.text = element_text(size = rel(1.25)),
        panel.grid.major.y = element_line(colour = "grey90"),
        panel.ontop = FALSE,
        text=element_text(size=18),
        axis.line = element_line(size = .5, colour = "black"))
ggsave("figures/russia-traj.pdf", width=8, height = 6, units = "in")

#####################################################################

#### Figure 2: Comparing headline target to change in emissions growth rate

dhsp <- read_dta("data/dhsp-v4-july.dta")

# Nominal to change in growth rate
ggplot(dhsp, aes(paris_nominal, delta_indcgrowthrate), label=iso) + 
  # Add smooth fit
  geom_smooth() + 
  # Add marker points at coordinates
  geom_text(aes(label=iso), hjust=0, vjust=0) +
  # Settings
  scale_x_continuous(breaks = seq(-20,60,10)) +
  scale_y_continuous(breaks = seq(-10,15,5), limits=c(-12,15)) +
  labs(x="Headline GHG reduction", y="Change in GHG growth rate under INDC") +
  annotate("text", x=-18, y=15, label = "Slower emissions growth under INDC", adj=0) +
  annotate("text", x=-18, y=-11, label = "Faster emissions growth under INDC", adj=0) +
  annotate("text", x=37, y=-12, label = "Larger headline reduction in INDC", adj=0) +
  annotate("text", x=0, y=-12, label = "Smaller headline reduction in INDC", adj=0) +
  theme(legend.position = "none",
        panel.background = element_rect(fill=NA), 
        axis.text = element_text(size = rel(1.25)),
        panel.grid.major.x = element_line(colour = "grey90"),
        panel.grid.major.y = element_line(colour = "grey90"),
        panel.ontop = FALSE,
        text=element_text(size=18),
        axis.line = element_line(size = .5, colour = "black"))
ggsave("figures/nominal-growthrate.pdf", width=8, height = 6, units = "in")


#####################################################################

#### Figure 3: Comparing headline target to range of possible percentage reductions
dhsp <- read_dta("data/dhsp-v4-july.dta")
dhsp$paris_baseyeartarget[dhsp$iso=="BEN"] <- 0 # Benin doesn't have a base year target

# Opportunism of the reference year, versus range of possible reference points
dhsp %>%
  # Re-shape data to calculate minimum and maximum values of delta_`yr'
  filter(!is.na(delta_2000)) %>%
  select(iso, paris_nominal, paris_baseyeartarget, starts_with("delta"), -delta_indcgrowthrate) %>%
  gather(key = key, value = value, starts_with("delta")) %>%
  group_by(iso) %>%
  # Calculate min and max
  mutate(max = round(max(value, na.rm=T), 2)) %>%
  mutate(min = round(min(value, na.rm=T), 2)) %>%
  ungroup() %>%
  spread(key, value) %>%
  select(iso, paris_nominal, paris_baseyeartarget, min, max) %>%
  # Create a rank variable to make plotting easier
  mutate(mid = (min+max)/2) %>%
  arrange(mid) %>%
  mutate(mid_rank = seq(1:nrow(.))) %>%
  # Re-label target types
  mutate(Target = ifelse(paris_baseyeartarget==1, "Base year", "BAU scenario")) %>%
  # -- Plot -- 
  ggplot(., aes(mid_rank, I(paris_nominal), colour=Target)) +
  # Set background line at zero
  geom_hline(yintercept = 0, linetype=2) +
  # Set marker points for base year and bau targets and plot
  scale_shape_manual(values = c(16, 17)) +
  geom_point(size=2.25, aes(shape=Target)) +
  # Add range bars and end points
  scale_color_manual(values = c("black", "blue")) +
  geom_segment(aes(x = mid_rank, xend = mid_rank, y=min, yend=max)) +
  geom_point(aes(mid_rank, min), pch="-", size=4.5) +
  geom_point(aes(mid_rank, max), pch="-", size=4.5) +
  # Settings
  coord_cartesian(ylim= c(-200, 100)) +
  scale_y_continuous(breaks=seq(-200,100,25)) +
  scale_x_continuous(breaks=NULL) +
  labs(y="Percentage emissions reduction", x="") +
  theme(legend.position = "bottom",
        panel.background = element_rect(fill=NA), 
        axis.text = element_text(size = rel(1.25)),
        panel.grid.major.y = element_line(colour = "grey90"),
        panel.ontop = FALSE,
        text=element_text(size=18),
        axis.line = element_line(size = .5, colour = "black"))
ggsave("figures/nominal-minmax.pdf", width=8, height = 6, units = "in")

#####################################################################

#### Figure 4: Coefficients from different delta_`yr' values

## Load in the data and re-shape

# Awareness variable
drummond_boot_awareness <- read.delim("figures/drummond_boot_awareness.txt", header=FALSE)
awareness <- drummond_boot_awareness[4:8,2:ncol(drummond_boot_awareness)]
rownames(awareness) <- c("beta", "se", "ci_low", "p_val", "ci_high")
awareness <- as.data.frame(t(awareness))
awareness <- mutate_all(awareness, function(x) as.numeric(as.character(x)))
awareness$year <- 1990:2014

# Serious variable
drummond_boot_serious <- read.delim("figures/drummond_boot_serious.txt", header=FALSE)
serious <- drummond_boot_serious[4:8,2:ncol(drummond_boot_serious)]
rownames(serious) <- c("beta", "se", "ci_low", "p_val", "ci_high")
serious <- as.data.frame(t(serious))
serious <- mutate_all(serious, function(x) as.numeric(as.character(x)))
serious$year <- 1990:2014

# Plot awareness
ggplot(awareness, aes(year, beta)) +
  geom_hline(yintercept = 0, linetype=2) +
  geom_segment(aes(x=year, xend=year, 
                   y=ci_low, yend=ci_high), size=1.15) +
  geom_point(size=2) +
  # geom_hline(yintercept = 0.237, col="red") + # Model II point estimate
  # geom_hline(yintercept = 0.011, col="red", linetype=2) + # Model II lower CI
  # geom_hline(yintercept = 0.462, col="red", linetype=2) + # Model II upper CI
  scale_y_continuous(limits=c(-2,4.3), breaks=(-2:4)) +
  scale_x_continuous(breaks=c(1990, 1995, 2000, 2005, 2010, 2014)) +
  labs(y="Coefficient and 95% confidence interval", x="Year", title="Awareness") +
  theme(legend.position = "none",
        panel.background = element_rect(fill=NA), 
        axis.text = element_text(size = rel(1.25)),
        panel.grid.major.x = element_line(colour = "grey90"),
        panel.grid.major.y = element_line(colour = "grey90"),
        panel.ontop = FALSE,
        text=element_text(size=18),
        axis.line = element_line(size = .5, colour = "black")) -> gg_awareness
gg_awareness

# Plot serious
ggplot(serious, aes(year, beta)) +
  geom_hline(yintercept = 0, linetype=2) +
  geom_segment(aes(x=year, xend=year, 
               y=ci_low, yend=ci_high), size=1.15) +
  geom_point(size=2) +
  # geom_hline(yintercept = 0.022, col="red") + # Model II point estimate
  # geom_hline(yintercept = -0.128, col="red", linetype=2) + # Model II lower CI
  # geom_hline(yintercept = 0.192, col="red", linetype=2) + # Model II upper CI
  scale_y_continuous(limits=c(-2,4.3), breaks=(-2:4)) +
  scale_x_continuous(breaks=c(1990, 1995, 2000, 2005, 2010, 2014)) +
  labs(y="", x="Year", title = "Perceived threat") +
  theme(legend.position = "none",
        panel.background = element_rect(fill=NA), 
        axis.text = element_text(size = rel(1.25)),
        panel.grid.major.x = element_line(colour = "grey90"),
        panel.grid.major.y = element_line(colour = "grey90"),
        panel.ontop = FALSE,
        text=element_text(size=18),
        # text=element_text(size=18, family = "LM Sans 10"), # LM Sans 10
        axis.line = element_line(size = .5, colour = "black")) -> gg_serious
gg_serious

## Combine plots and save
coef_plot <- gridExtra::grid.arrange(gg_awareness, gg_serious, nrow = 1, ncol=2)
ggsave("figures/years-coefplot.pdf", coef_plot, width = 14, height = 6, units = "in")

#####################################################################

#### Figure 5: Favorable base years

## Load in the data and re-shape
dhsp <- read_dta("data/dhsp-v4-july.dta")

dhsp %>%
  filter(paris_baseyeartarget==1, 
         !is.na(paris_nominal)) %>%
  select(iso, paris_nominal, paris_baseyeartarget, starts_with("delta"), -delta_indcgrowthrate) %>%
  select(-iso, -paris_baseyeartarget) %>%
  t(.) -> df
# Countries as columns, delta_`yr' as rows and paris_nominal (paris_headline) as first entry

rank <- rep(NA, ncol(df))
for (i in 1:length(rank)){ 
  rank[i] <- rank(df[,i], ties.method = "min")
}  
rank <- (-1*rank)+27 # So that most favourable years have low values

ggplot(as.data.frame(rank), aes(rank)) + 
  geom_histogram(breaks=seq(0, 26, by=1), color="white", fill="grey30") + 
  scale_x_continuous(breaks=seq(0,25, 5)) +
  scale_y_continuous(breaks=seq(0,14,2)) +
  labs(y="Count of countries", 
       x = "Favourability rank of a country's chosen base year \nCountries with base year targets only") +
  theme(legend.position = "none",
        panel.background = element_rect(fill=NA), 
        axis.text = element_text(size = rel(1.25)),
        panel.grid.major.x = element_line(colour = "grey90"),
        panel.grid.major.y = element_line(colour = "grey90"),
        panel.ontop = FALSE,
        text=element_text(size=18),
        axis.line = element_line(size = .5, colour = "black"))
ggsave("figures/rank.pdf", width=8, height = 6, units = "in")


#####################################################################

#### Table #: Common samples ####

## Load in and reshape output tables individually

common123 <- read.delim("tables/dhsp_commonsample123.txt", header=FALSE)
common123_aware <- common123[c(4:8), c(2,4,6)]
colnames(common123_aware) <- c("dhsp", "nominal", "delta")
rownames(common123_aware) <- c("beta", "se", "ci_low", "p_val", "ci_high")
common123_aware <- as.data.frame(t(common123_aware))
common123_aware$sample <- "1, 2, 3"
common123_aware$outcome <- row.names(common123_aware)

common123_serious <- common123[c(9:13), c(3,5,7)]
colnames(common123_serious) <- c("dhsp", "nominal", "delta")
rownames(common123_serious) <- c("beta", "se", "ci_low", "p_val", "ci_high")
common123_serious <- as.data.frame(t(common123_serious))
common123_serious$sample <- "1, 2, 3"
common123_serious$outcome <- row.names(common123_serious)

common23 <- read.delim("tables/dhsp_commonsample23.txt", header=FALSE)
common23_aware <- common23[c(4:8), c(2,4,6)]
colnames(common23_aware) <- c("dhsp", "nominal", "delta")
rownames(common23_aware) <- c("beta", "se", "ci_low", "p_val", "ci_high")
common23_aware <- as.data.frame(t(common23_aware))
common23_aware$sample <- "2, 3"
common23_aware$outcome <- row.names(common23_aware)

common23_serious <- common23[c(9:13), c(3,5,7)]
colnames(common23_serious) <- c("dhsp", "nominal", "delta")
rownames(common23_serious) <- c("beta", "se", "ci_low", "p_val", "ci_high")
common23_serious <- as.data.frame(t(common23_serious))
common23_serious$sample <- "2, 3"
common23_serious$outcome <- row.names(common23_serious)

## Combine intermediate tables into one summary regression table

common_aware <- rbind(common123_aware, common23_aware)
common_serious <- rbind(common123_serious, common23_serious)
common <- cbind(common_aware, common_serious)  
outputs <- common[,c(6, 7, 1, 2, 4, 8, 9, 11)]  
outputs <- dplyr::arrange(outputs, outcome, sample)  
# Set number of observations
outputs$obs <- ifelse(outputs$sample=="1, 2, 3", 35,
                      ifelse(outputs$sample=="2, 3" & outputs$outcome!="dhsp", 60, 35))
# Latex table output
stargazer(outputs, summary = F, rownames = F)


#####################################################################

#### Table #: Summary of key variables for all states ####

## Load in and reshape data

dhsp <- read_dta("data/dhsp-v4-july.dta")
dhsp$paris_baseyeartarget[dhsp$iso=="BEN"] <- 0 # Benin doesn't have a base year target

dhsp %>%
  select(iso, paris_ghgtarget, paris_nonghgtarget, paris_actions,
         paris_baseyeartarget, paris_scenariotarget, paris_intensitytarget, paris_trajectorytarget, paris_fixedleveltarget,
         dhsp = unconditional, headline = paris_nominal, delta_indcgrowthrate, paris_referenceyear, lee_awareness,
         model1, model2, model3, model4) %>%
  mutate(target = ifelse(paris_nonghgtarget==1 | paris_actions==1, "Non-GHG target", 
                         ifelse(paris_trajectorytarget == 1 | 
                                  paris_fixedleveltarget == 1 | 
                                  paris_intensitytarget == 1, "Other GHG target",
                         ifelse(paris_scenariotarget ==1, "BAU target", 
                         ifelse(paris_baseyeartarget == 1, "Base year target", NA))))) %>%
  mutate(models = ifelse(model1==1 & model2==1 & model3==1, "1, 2, 3", 
                         ifelse(model1==1 & model2==1 & model3 == 0, "1, 2",
                         ifelse(model1==1 & model2==0 & model3 == 0, "1",
                         ifelse(model1==0 & model2==1 & model3 == 1, "2, 3",
                         ifelse(model1==0 & model2==1 & model3 == 0, "2",
                         ifelse(model1==0 & model2==0 & model3 == 0, ".", NA))))))) %>%
  mutate(dhsp = round(dhsp, 1)) %>%
  mutate(headline = round(headline, 1)) %>%
  mutate(delta_indcgrowthrate = round(delta_indcgrowthrate, 1)) %>% 
  mutate(delta_indcgrowthrate = round(delta_indcgrowthrate, 1)) %>% 
  mutate(lee_awareness = round(lee_awareness, 1)) %>%
  arrange(iso) %>%
  select(iso, target, dhsp, headline, delta_indcgrowthrate, 
         baseyear = paris_referenceyear, awareness = lee_awareness, models) -> key_variables
key_variables[is.na(key_variables)] <- "."

stargazer(key_variables, summary = F, rownames = F)





